function or_column_map

global P net

samples = 1000;
pixels = 10;



P = zeros(pixels^2,samples);
size(P)
for i = 1:samples
    cycles = 2;
    xphase = (-1+2*rand)*cycles*pi;
    yphase = (-1+2*rand)*cycles*pi;
    angle = rand*pi;
    contrast = 1;
    sd_x = 1;
    sd_y=4;

    g=contrast * gaborfilter(pixels,sd_x*pi,sd_y*pi,angle,cycles,xphase,yphase);

    if mod(i,100)==0
        imagesc(g,[-1 1])
        colormap(gray);
        drawnow;
    end

    size(P)
    size(g(:))
    P(:,i)=g(:);


    %pause(.5);    
    
end



% Now train a SOM on the data samples
net = newsom(P,[10 10]);  % this is a 2D topography network.

%%
% We can visualize the network we have just created with PLOTSOM.
% 
% Each neuron is represented by a red dot at the location of its two weights.
% Initially all the neurons have the same weights in the middle of the vectors,
% so only one dot appears.

plotsom(net.iw{1,1},net.layers{1}.distances)

%%
% Now we train the map on the 1000 vectors for 1 epoch and replot the network
% weights.
% 
% After training, note that the layer of neurons has begun to self-organize so
% that each neuron now classifies a different region of the input space, and
% adjacent (connected) neurons respond to adjacent regions.

net.trainParam.epochs = 100;
net = train(net,P);

w=net.iw{1,1};
%plotsom(net.iw{1,1},net.layers{1}.distances)

figure(2)
k=1;
for i = 1:10
    for j = 1:10
        subplot(10,10,k);
        imagesc(reshape( w(k,:),pixels,pixels ),[-1 1]);
        colormap(gray);
        k=k+1;
    end
end


%% 
% Now test the network









% Modified by Jeff McKinstry, 10-28-09
% Allows passing in an array of x values for the cross-section of the Gabor
% function.  This allows you to set the phase and the number of pixels in
% your output array.  
% 
% Sample usage:
% imagesc(gaborfilter1(-2*pi:pi/4:2*pi,pi,10*pi,pi/4)),colormap(gray)
%
%The Gabor filter is basically a Gaussian (with variances sx and sy along x and y-axes respectively)
%modulated by a sinusoid (with centre frequencies U and V along x and y-axes respectively) 
%described by the following equation
%%
%                            -1     x' ^     y'  ^             
%%% G(x,y,theta,f) =  exp ([----{(----) 2+(----) 2}])*cos(2*pi*f*x');
%                             2    sx'       sy'
%%% x' = x*cos(theta)+y*sin(theta);
%%% y' = y*cos(theta)-x*sin(theta);

%% Describtion :

%% one_d_fun : Input array of x values for the cross-section of the Gabor function; if n==length(one_d_fun) then G is nXn.
%% Sx & Sy : Variances along x and y-axes respectively
%% theta : The orientation of Gabor filter

%% G : The output filter as described above



%%  Author : Ahmad poursaberi  e-mail : a.poursaberi@ece.ut.ac.ir
%%          Faulty of Engineering, Electrical&Computer Department,Tehran
%%          University,Iran,June 2004

function G = gaborfilter(pixels,Sx,Sy,theta,cycles,xphase,yphase);

step = 2*pi/((pixels-1)/cycles);
one_d_xfun = (-cycles*pi+xphase):step:(cycles*pi+xphase);
one_d_yfun = (-cycles*pi+yphase):step:(cycles*pi+yphase);


% rotate
[x y] = meshgrid(one_d_xfun,one_d_yfun);
xPrime = x * cos(theta) + y * sin(theta);
%imagesc(xPrime)
yPrime = y * cos(theta) - x * sin(theta);
%figure;imagesc(yPrime)
G = exp(-.5*((xPrime/Sx).^2+(yPrime/Sy).^2)).*cos(xPrime);

%G=double(real(G));

